# Set up ---------------------------------
# Load Libraries 
# library(groundhog)
# groundhog.library(tidyverse, "2023-01-01")
# groundhog.library(lfe, "2023-01-01")
# groundhog.library(fixest, "2023-01-01")
# groundhog.library(egg, "2023-01-01")
# groundhog.library(scales, "2023-01-01")
library(tidyverse)
library(lfe)
library(fixest)
library(egg)
library(scales)

adl <- read_csv("output/data/adl_full_subset.csv")
adl2013 <- read_csv("output/data/adl2013_outgroup_attitudes.csv")
country_cov <- read_csv("input/predictors-of-support.csv")
load("input/UN_signatories.RData")


# Functions -----------------
glance.felm_custom <- function(x, ...) {
  out <- broom:::glance.felm(x, ...)
  for (fe in names(x$fe)) {
    out[[paste('Fixed effects:', fe)]] <- 'X'
  }
  out
}
felm_custom <- function(...) {
  out <- felm(...)
  class(out) <- c('felm_custom', class(out))
  out
}

run_baseline_mod <- function(dv, df) {
  dv <- as.symbol(dv)
  eval(
    bquote(
      felm_custom(
        .(dv) ~ education*signed_num | wave + country | 0 | country,
        data = df
      )
    )
  )
}

run_control_mod <- function(dv, df) {
  dv <- as.symbol(dv)
  eval(
    bquote(
      felm_custom(
        .(dv) ~ education*signed_num + economy + persfinance + 
          polsit + gender + yearborn + newssource | 
          wave + country | 0 | country,
        data = df
      )
    )
  )
}

print_paran <- function(text) {
  paste0("(", text, ")")
}

print_star_coef <- function(mod, var) {
  df <- tidy(mod)
  var_row <- df[df$term == var,]
  coef <- round(var_row$estimate, 3)
  stars <- 
    if (var_row$p.value < 0.001) {
      "***"
    } else if (var_row$p.value < 0.01) {
      "**"
    } else if (var_row$p.value < 0.05) {
      "*"
    } else ""
  paste0(coef, stars)
}


# Preparing data for analysis ---------

# Change outgroup sentiment and stereotypes to linear categories
# Change evaluation of sentiment from in-groups to be NA (i.e. Christians don't rate other Christians)
adl_num <- 
  adl %>% 
  select(-tfjneighbor) %>% # inconsistent data
  mutate(
    across(ends_with("spop"), recode,
           `favorable` = 0, 
           `recognize, can't rate` = 1,
           `unfavorable` = 2)
  ) %>% 
  mutate(
    jewspop = case_when(
      religion == "jewish" ~ NA_real_,
      TRUE ~ jewspop
    ),
    muslimspop = case_when(
      religion == "muslim" ~ NA_real_,
      TRUE ~ muslimspop
    ),
    christianspop = case_when(
      religion == "christian" ~ NA_real_,
      TRUE ~ christianspop
    ),
    hinduspop = case_when(
      religion == "hindu" ~ NA_real_,
      TRUE ~ hinduspop
    ),
    buddhistspop = case_when(
      religion == "buddhist" ~ NA_real_,
      TRUE ~ buddhistspop
    ),
    wave_region = paste0(wave, ":", world_region),
    signed_num = factor(signed_num, levels = c(0,1,2))
  ) 

adl_streotypes <- 
  adl %>% 
  select(-tfjneighbor) %>% 
  rowwise() %>% 
  mutate(
    count_t = 
      sum(
        c_across(starts_with("tf")) == "probably true", 
        na.rm = T
      ),
    count_f = 
      sum(
        c_across(starts_with("tf")) == "probably false", 
        na.rm = T
      ),
    count_dk = 
      sum(
        c_across(starts_with("tf")) == "don't know", 
        na.rm = T
      ),
    net_t = count_t - count_f,
    prop_t = count_t / 11,
    prop_f = count_f / 11,
    prop_dk = count_dk / 11,
    prop_net_t = net_t / 11
  ) %>% 
  ungroup()

adl_streotypes <- 
  adl_streotypes %>% 
  mutate(
    wave_region = paste0(wave, ":", world_region),
    signed_num = factor(signed_num, levels = c(0,1,2))
  )

# Clean just 2013 data 
un_countries <- unique(adl2013$country)
un_both <- intersect(un07, un15)
un_either <- setdiff(union(un07, un15), un_both)
adl2013 <- 
  adl2013 %>% 
  mutate(
    signed_un = case_when(
      country %in% un_both ~ "Supported both",
      country %in% un_either ~ "Supported Holocaust only",
      TRUE ~ "Supported none"
    ) %>% 
      factor(levels = c("Supported none", "Supported Holocaust only", "Supported both"))
  ) %>% 
  rowwise() %>%
  mutate(
    count_t = 
      sum(
        c_across(starts_with("tf")) == "probably true", 
        na.rm = T
      )
  ) %>% 
  ungroup()

# Change variables names and set levels 
adl2013 <- 
  adl2013 %>% 
  mutate(
    ED_Age = recode(ED_Age, `[vol] refused` = "DK/refused"),
    persfinance = recode(persfinance, `don't knowrefused` = "DK/refused"),
    economy = recode(economy, `don't know/refused` = "DK/refused"),    
    polsit = recode(polsit, `don't know/refused` = "DK/refused"),    
    newssource = recode(newssource, `don't know/refused` = "DK/refused"),   
    relattend = recode(relattend, `don't know/refused` = "DK/refused"),
    religion = recode(religion, 
                      jewish = "Jewish", christian = "Christian", 
                      muslim = "Muslim", buddhist = "Buddhist",
                      hindu = "Hindu", 
                      `none/atheist/agnostic` = "None/Atheist/Agnostic",
                      `don't know/refused` = "DK/refused"),
    ED_Age = factor(ED_Age, 
                    levels = c("5-12", "13-18", "19-22", 
                               "23+", "DK/refused")),
    religion = factor(religion, 
                      levels = c("Jewish", "Christian", "Muslim", 
                                 "Buddhist", "Hindu", "None/Atheist/Agnostic",
                                 "DK/refused")),
    relattend = factor(relattend,
                       levels = c("less than once a month", "one to two times a month",
                                  "once a week", "multiple times a week", 
                                  "daily", "DK/refused")),
    persfinance = factor(persfinance, 
                         levels = c("poor", "not so good",
                                    "good", "excellent", "DK/refused")),
    economy = factor(economy, 
                     levels = c("poor", "not so good",
                                "good", "excellent", "DK/refused")),
    polsit = factor(polsit, 
                    levels = c("very unstable", "somewhat unstable",
                               "somewhat stable", "very stable", "DK/refused")),
    newssource = factor(newssource, 
                        levels = c("other people", "newspapers", "radio",
                                   "television", "the internet", "DK/refused"))
  )


# Country covariates 
# Merge region indicators 
cntry_region <- 
  adl2013 %>% 
  distinct(country, world_region, dominant_relig, hdi)
country_cov <- 
  country_cov %>% 
  left_join(cntry_region, by = "country")

# Create variables 
country_cov <- 
  country_cov %>% 
  mutate(
    sign_1 = if_else(signed_07 == 1 | signed_15 == 1, 1L, 0L),
    sign_2 = if_else(signed_07 == 1 & signed_15 == 1, 1L, 0L)
  ) 

# Rename long variables 
names(country_cov)[7:11] <- 
  c("jew_pop", "regime", "gdp", "islam_death_79_00", "islam_death_01_12")

# Fix variables
country_cov <- 
  country_cov %>% 
  mutate(
    jew_pop = as.numeric(jew_pop),
    regime = as.numeric(regime),
    gdp = as.numeric(gdp),
    dominant_relig = factor(dominant_relig, 
                            levels = c("other", "christian", "muslim", 
                                       "buddhist", "hindu")),
    world_region = factor(world_region,
                          levels = c("oceania", "americas", "e europe",
                                     "w europe", "mena", "sub-saharan africa",
                                     "asia")),
    hdi = factor(hdi, levels = c("1-47", "142-187", "48-141"))
  )


# Analysis ----------------------------

# Figure 1 
# Jewish stereotypes 
# plot means for wave 1
count_by_edu <- 
  adl_streotypes %>% 
  ungroup() %>% 
  filter(wave == 1) %>% 
  mutate(
    sum_t = prop_t * 11,
    sum_f = prop_f * 11,
    sum_dk = prop_dk * 11,
    signed_num = recode(
      signed_num, `0` = "Supported none", 
      `1` = "Supported Holocaust only", 
      `2` = "Supported both")
  ) %>% 
  group_by(education, signed_num) %>% 
  summarise(
    mean_count_t = mean(sum_t, na.rm = T),
    mean_count_f = mean(sum_f, na.rm = T),
    mean_count_dk = mean(sum_dk, na.rm = T)
  ) %>% 
  filter(!is.na(education))

edu_means_plot <- 
  count_by_edu %>% 
  pivot_longer(
    cols = contains("_count_"),
    names_to = "who",
    values_to = "count"
  ) %>%
  mutate(
    prop = 100 * count / 11,
    who = factor(who, levels = c("mean_count_t", "mean_count_f", "mean_count_dk"))
  ) %>% 
  ggplot(
    aes(x = education, y = prop)
  ) +
  theme_bw() + 
  geom_point(aes(color = who, shape = who)) + 
  geom_line(aes(color = who, linetype = who)) +
  ggtitle("") + 
  scale_x_continuous(
    name = "Age at final year of education",
    breaks = c(1, 2, 3, 4),
    labels = c(expression(""<=12), "13-18", "19-22", "23+") 
  ) +
  scale_y_continuous(
    name = "Mean responses",
    labels = function(x) paste0(x, "%"),
    limits = c(0, 100)
  ) +
  scale_color_discrete(
    name = NULL,
    labels = c("Probably true", "Probably false", "Don't know/recognize")
  ) +
  scale_linetype_discrete(
    name = NULL,
    labels = c("Probably true", "Probably false", "Don't know/recognize")
  ) +
  scale_shape_discrete(
    name = NULL,
    labels = c("Probably true", "Probably false", "Don't know/recognize")
  ) +
  facet_grid(cols = vars(signed_num))


# Table 1a (favorability regressions)

relig_favor <- 
  tibble(
    title = c("Jews", "Christians", "Muslims", 
              "Buddhists", "Hindus"),
    who = c("jewish", "christian", "muslim", 
            "buddhist", "hindu"),
    variable = c("jewspop", "christianspop", "muslimspop", 
                 "buddhistspop", "hinduspop")
  )
relig_control_mods <- vector("list", length = nrow(relig_favor))
for(i in 1:nrow(relig_favor)) {
  df <- filter(adl_num, religion != relig_favor$who[i])
  df$signed_num <- factor(df$signed_num, levels = c(0,1,2))
  relig_control_mods[[i]] <- 
    run_control_mod(relig_favor$variable[i], df)
}

# Loop through to change baseline for marginal effects 
fav_ME_mods <- vector("list", length = 2*nrow(relig_favor))
for(i in 1:nrow(relig_favor)) {
  df <- filter(adl_num, religion != relig_favor$who[i])
  df$signed_num <- factor(df$signed_num, levels = c(1,2,0))
  fav_ME_mods[[2*i - 1]] <- 
    run_control_mod(relig_favor$variable[i], df)
  df$signed_num <- factor(df$signed_num, levels = c(2,1,0))
  fav_ME_mods[[2*i]] <- 
    run_control_mod(relig_favor$variable[i], df)
}


# Table 1b (stereotype regressions)
stereotype_vars <- 
  tibble(
    title = c("Pct true", "Pct false", "Net endorse"),
    variable = c("prop_t", "prop_f", "prop_net_t")
  )
stereotype_control_mods <- vector("list", length = nrow(stereotype_vars))
for(i in 1:nrow(stereotype_vars)) {
  df$signed_num <- factor(df$signed_num, levels = c(0,1,2))
  stereotype_control_mods[[i]] <-
    run_control_mod(stereotype_vars$variable[i], adl_streotypes)
}

stereotype_ME_mods <- vector("list", length = 2*nrow(stereotype_vars))
for(i in 1:nrow(stereotype_vars)) {
  df <- adl_streotypes
  df$signed_num <- factor(df$signed_num, levels = c(1,2,0))
  stereotype_ME_mods[[2*i - 1]] <- 
    run_control_mod(stereotype_vars$variable[i], df)
  df$signed_num <- factor(df$signed_num, levels = c(2,0,1))
  stereotype_ME_mods[[2*i]] <- 
    run_control_mod(stereotype_vars$variable[i], df)
}


# Figure A1 & A2 
vars <- 
  c("prop_t", "prop_f", "jewspop", "muslimspop", 
    "christianspop", "hinduspop", "buddhistspop")
var_title <- 
  c("Proportion antisemitic stereotype true", 
    "Proportion antisemitic stereotype false", 
    "Favorability towards Jews", 
    "Favorability towards Muslims", 
    "Favorability towards Christians", 
    "Favorability towards Hindus", 
    "Favorability towards Buddhists")
outcome_hist <- vector("list", length = length(vars))

adl_streotypes <- mutate(adl_streotypes, wave_chr = paste0("Wave ", wave))

for(v in 1:length(vars)) {
  
  if(v <= 2) {
    
    outcome_hist[[v]] <- 
      adl_streotypes %>%
      ggplot(aes(x = !!sym(vars[v]))) + 
      theme_bw() + 
      geom_histogram(
        binwidth = 0.1
      ) +
      xlab(NULL) + ylab(NULL) +
      ggtitle(var_title[v]) + 
      facet_wrap(vars(wave_chr), nrow = 1)
    
  } else {
    outcome_hist[[v]] <- 
      adl_streotypes %>%
      filter(!is.na(!!sym(vars[v]))) %>% 
      ggplot(aes(x = !!sym(vars[v]))) + 
      theme_bw() + 
      geom_bar() +
      xlab(NULL) + ylab(NULL) +
      ggtitle(var_title[v]) + 
      scale_x_discrete(labels= c("Favorable", "Can't rate", "Unfavorable")) +
      facet_wrap(vars(wave_chr), nrow = 1)
  }
  
  
}

outcome_hist_pct <- vector("list", length = length(vars))

for(v in 1:length(vars)) {
  
  if(v <= 2) {
    
    outcome_hist_pct[[v]] <- 
      adl_streotypes %>%
      ggplot(aes(x = !!sym(vars[v]))) + 
      theme_bw() + 
      geom_histogram(
        aes(y = after_stat(density)*0.1), 
        binwidth=0.1
      ) +
      xlab(NULL) + ylab(NULL) +
      ggtitle(var_title[v]) + 
      facet_wrap(vars(wave_chr), nrow = 1)
    
  } else {
    outcome_hist_pct[[v]] <- 
      adl_streotypes %>%
      filter(!is.na(!!sym(vars[v]))) %>% 
      group_by(wave_chr, !!sym(vars[v])) %>% 
      summarise(n = n()) %>% 
      group_by(wave_chr) %>% 
      mutate(
        wave_total = sum(n),
        pct = n / wave_total
      ) %>% 
      ungroup() %>% 
      ggplot(aes(x = !!sym(vars[v]), y = pct)) + 
      theme_bw() + 
      geom_bar(stat="identity") +
      xlab(NULL) + ylab(NULL) +
      ggtitle(var_title[v]) + 
      scale_x_discrete(labels= c("Favorable", "Can't rate", "Unfavorable")) +
      facet_wrap(vars(wave_chr), nrow = 1)
  }
  
  
}

# A1 
stereo_dist <- 
  ggarrange(outcome_hist_pct[[1]], outcome_hist_pct[[2]], ncol = 1)

# A2
outgroup_plot <- 
  ggarrange(outcome_hist_pct[[3]], outcome_hist_pct[[4]], outcome_hist_pct[[5]],
            outcome_hist_pct[[6]], outcome_hist_pct[[7]], ncol = 1)


# Figure A3

vars <- c("count_t", "ED_Age", "gender", "age", "dominant_relig", "relattend")
var_lab <- c("N antisemitic stereotype true", 
             "Last education age", "Gender",
             "Age", "Dominant religion", "Religiosity")

# Add line breaks to variable names 
levels(adl2013$religion) <- 
  c("Jewish", "Christian", "Muslim", "Buddhist", "Hindu",
    "None\nAtheist\nAgnostic", "DK/refused")
levels(adl2013$relattend) <- 
  c("less than\nonce\na month", "1-2\ntimes\na month", "once\na week", 
    "multiple\ntimes\n a week", "daily",
    "DK/refused")
plot_rows <- vector("list", length = length(vars))

for(i in 1:length(vars)) {
  
  # Tally proportions
  pct_var <- 
    adl2013 %>% 
    group_by(signed_un) %>% 
    count(!!sym(vars[i])) %>% 
    mutate(
      prop = 100 * n / sum(n),
      var_name = var_lab[i]
    ) %>% 
    rename(category = !!sym(vars[i]))
  
  # Plot 
  plot_rows[[i]]  <- 
    pct_var %>% 
    filter(!is.na(category)) %>%
    ggplot(aes(x = category, y = prop)) + 
    theme_bw() + 
    geom_bar(stat="identity") + 
    xlab(NULL) + ylab(NULL) + 
    ggtitle(var_lab[i]) +
    scale_y_continuous(labels = label_percent(accuracy = 1, scale = 1)) +
    facet_grid(cols = vars(signed_un))
}

# A3 
compare_signage_vars <- 
  ggarrange(plots = plot_rows, ncol = 1)


# Figure A4(a) and A4(b)
adl2013 <- 
  adl2013 %>%
  mutate(
    probably_BA = if_else(ED_Age %in% c("19-22", "23+"), 
                          "Last education age 19+", "Last education age <=18", NA_character_)
  )

vars <- c("count_t", "gender", "age", "religion", "relattend")
var_lab <- c("N antisemitic stereotype true", 
             "Gender", "Age", "Religion", "Religiosity")

for(i in 1:length(vars)) {
  
  # Tally proportions
  pct_var <- 
    adl2013 %>% 
    group_by(signed_un, probably_BA) %>% 
    count(!!sym(vars[i])) %>% 
    mutate(
      prop = 100 * n / sum(n),
      var_name = var_lab[i]
    ) %>% 
    rename(category = !!sym(vars[i]))
  
  # Plot 
  plot_rows[[i]]  <- 
    pct_var %>% 
    filter(!is.na(category)) %>%
    ggplot(aes(x = category, y = prop)) + 
    theme_bw() + 
    geom_bar(stat="identity") + 
    xlab(NULL) + ylab(NULL) + 
    ggtitle(var_lab[i]) +
    scale_y_continuous(labels = label_percent(accuracy = 1, scale = 1)) +
    facet_grid(cols = vars(signed_un), rows = vars(probably_BA),
               switch = "y")
}

plot_a4a <- 
  ggarrange(plots = plot_rows[1:3], ncol = 1)
plot_a4b <- 
  ggarrange(plots = plot_rows[4:5], ncol = 1)


# Table A3 

# Jewish sentiment 
americas_jewspop <- 
  felm(
    jewspop ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_num, religion != "jewish" & world_region == "americas")
  ) 
w_europe_jewspop <- 
  felm(jewspop ~ education*signed_num + economy + persfinance + 
         polsit + gender + yearborn + newssource|wave + country |
         0 | country,
       data = filter(adl_num, religion != "jewish" & world_region == "w europe"))

e_europe_jewspop <- 
  felm(
    jewspop ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_num, religion != "jewish" & world_region == "e europe")
  )
oceania_jewspop <- 
  felm(
    jewspop ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_num, religion != "jewish" & world_region == "oceania")
  )
asia_jewspop <- 
  felm(
    jewspop ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_num, religion != "jewish" & world_region == "asia")
  )
mena_jewspop <- 
  felm(
    jewspop ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_num, religion != "jewish" & world_region == "mena")
  )
sub_africa_jewspop <- 
  felm(
    jewspop ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_num, religion != "jewish" & world_region == "sub-saharan africa")
  )
jewspop_interact <- 
  felm(
    jewspop ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave_region + country |
      0 | country,
    data = filter(adl_num, religion != "jewish")
  )

# Loop through to get marginal effects 
region <- c("americas", "w europe", "e europe", "oceania", "asia", 
            "mena", "sub-saharan africa", "")
fav_ME_mods <- vector("list", length = 2*length(region))

for(i in 1:length(region)) {
  
  if(i != 8) {
    
    df <- filter(adl_num, religion != "jewish" & world_region == region[i])
    
  } else df <- filter(adl_num, religion != "jewish")
  
  df$signed_num <- factor(df$signed_num, levels = c(1,2,0))
  fav_ME_mods[[2*i - 1]] <- 
    run_control_mod("jewspop", df)
  df$signed_num <- factor(df$signed_num, levels = c(2,1,0))
  fav_ME_mods[[2*i]] <- 
    run_control_mod("jewspop", df)
  
}


# Table A4
americas <- 
  felm(
    prop_net_t ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_streotypes, world_region == "americas")
  ) 
w_europe <- 
  felm(prop_net_t ~ education*signed_num + economy + persfinance + 
         polsit + gender + yearborn + newssource|wave + country |
         0 | country,
       data = filter(adl_streotypes, world_region == "w europe"))

e_europe <- 
  felm(
    prop_net_t ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_streotypes, world_region == "e europe")
  )
oceania <- 
  felm(
    prop_net_t ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_streotypes, world_region == "oceania")
  )
asia <- 
  felm(
    prop_net_t ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_streotypes, world_region == "asia")
  )
mena <- 
  felm(
    prop_net_t ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_streotypes, world_region == "mena")
  )
sub_africa <- 
  felm(
    prop_net_t ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave + country |
      0 | country,
    data = filter(adl_streotypes, world_region == "sub-saharan africa")
  )
stereotype_interact <- 
  feols(
    prop_net_t ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource | country + wave^world_region,
    data = adl_streotypes
  )
stereotype_interact <- 
  felm(
    prop_net_t ~ education*signed_num + economy + persfinance + 
      polsit + gender + yearborn + newssource| wave_region + country |
      0 | country,
    data = adl_streotypes
  )

stereo_ME_mods <- vector("list", length = 2*length(region))

for(i in 1:length(region)) {
  
  if(i != 8) {
    
    df <- filter(adl_streotypes, world_region == region[i])
    
  } else df <- adl_streotypes
  
  df$signed_num <- factor(df$signed_num, levels = c(1,2,0))
  stereo_ME_mods[[2*i - 1]] <- 
    run_control_mod("prop_net_t", df)
  df$signed_num <- factor(df$signed_num, levels = c(2,1,0))
  stereo_ME_mods[[2*i]] <- 
    run_control_mod("prop_net_t", df)
  
}

# Table A5 
signed_1 <- 
  lm(
    sign_1 ~ jew_pop + regime + log(gdp) + 
      hdi + dominant_relig + world_region, 
    data = country_cov
  )

signed_2 <- 
  lm(
    sign_2 ~ jew_pop + regime + log(gdp) + 
      hdi + dominant_relig + world_region, 
    data = country_cov
  )




